function [parcube,b,timepoints,nii] = fmripod_parcube(theta,hdrint,hdrdbl)
% fmripod_parcube
%
%    [parcube,b]                = fmripod_parcube(theta,hdrint)
%    [parcube,b,timepoints,nii] = fmripod_parcube(theta,hdrint,hdrdbl)
%
% input:
%
%         theta : the parameters that come out of FMRIPOD if it was run with fold=1
%        hdrint : comes out of FMRIPOD.
%        hdrdbl : comes out of FMRIPOD.
%
% output:
%
%       parcube : Pseudo parameter array (Dx,Dy,Dz,Dt)
%             b : Logistic argument offset (see [1]).
%    timepoints : Temporal samplings  (requires hdrdbl given as input)
%           nii : Volume for nifti toolbox visualization (requires hdrdbl given as input)
%
% Copyright: Mads Dyrholm
b = theta(end);
theta = theta(1:end-1);
% hdrint fetch
Dx = hdrint(1);
Dy = hdrint(2);
Dz = hdrint(3);
fold = hdrint(12);
if fold>1
  error('You need to run fmripod with fold=1');
end
mode = hdrint(15)
R = hdrint(16);
Dt = hdrint(17);
% do it
if mode==1
  parcube = reshape(theta,[Dx,Dy,Dz,Dt]);
elseif mode==2
  parcube = zeros(Dx,Dy,Dz,Dt);
  rridx = 0
  for r = 1:R
    ux = theta(rridx+(1:Dx*Dy*Dz));
    for t = 1:Dt
      parcube(:,:,:,t) = parcube(:,:,:,t) + reshape(ux,[Dx,Dy,Dz])  * theta(rridx+Dx*Dy*Dz+t);
    end
    rridx = rridx + Dx*Dy*Dz + Dt;
  end
elseif mode==3
  parcube = zeros(Dx,Dy,Dz);
  rridx = 0
  for r = 1:R
    ux = theta(rridx+(1:Dx));
    uy = theta(rridx+Dx+(1:Dy));
    uz = theta(rridx+Dx+Dy+(1:Dz));
    for z = 1:Dz
      parcube(:,:,z) = parcube(:,:,z) + ux*uy'*uz(z);
    end
    rridx = rridx + Dx+Dy+Dz;
  end
end
% nii
if nargout>2
  timepoints = linspace(hdrdbl(14),hdrdbl(14)+hdrdbl(15)*(hdrint(17)-1),hdrint(17));
end
if nargout>3
  nii = make_nii(parcube, [hdrdbl([2 3 4])]);
end
